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1.      Introduction 


An  interest  in  the  properties  of  incompressible  fluid  flow  in  pipes  can  undoubtedly  be 
traced  as  far  back  as  the  era  of  the  Roman  empire.  A  brief  synopsis  of  the  historical 
development  of  fluid  mechanics  can  be  found  in  many  of  the  introductory  level  texts 
[Vennard  (1961)  for  example]. 

The  problem  of  accurately  solving  for  the  parameters  in  incompressible  fluid  flow 
in  piping  networks  is  dependent  upon  the  single-pipe  model  for  fluid  flow.  Many 
semi- empirical  methods  for  characterizing  the  flow  properties  have  appeared  in  the 
literature  [Matthew  (1981),  Wuori  (1993),  and  Churchill  (1977)].  Since  the  early  work 
by  Feigenbaum  (1978),  see  for  example  [Hofstadter  (1981)],  on  chaos  in  simple  sys- 
tems, considerable  progress  has  been  made  on  analytically  characterizing  turbulence 
in  fluids  [Landau  and  Lifshitz  (1987)].  Nevertheless,  many  engineering  problems  are 
solved  using  the  Stanton  diagram,  which  shows  Nikuradse's  measured  data. 

For  engineering  problems  which  require  repetitive  numerical  evaluation,  an  empir- 
ical relationship  facilitates  the  task  and  various  investigators  [Tsang  and  Kee  (1987), 
Round  (1980,  1985),  and  Colebrook  (1939)]  have  proposed  such  formulae.  These  for- 
mulae have  the  drawback  that  the  transition  region  between  laminar  and  turbulent 
flow  is  inadequately  represented.  Only  the  fluid  flow  properties  in  the  laminar  region 
are  fairly  well  understood  [Vennard  (1961)]  using  methods  of  fluid  analysis. 

One  of  the  earliest  and  most  well  known  methods  for  solving  pipeline  networks 
is  attributed  to  Hardy  Cross  [Giles  (1962)  and  Lindeburg  (1982)].  The  Hardy  Cross 
scheme,  not  only  requires  an  a  priori  guess  on  the  initial  flows,  but  upon  comparison 
with  more  recently  reported  methods  [Carnahan  and  Christensen  (1972),  Bending 
and  Hutchinson  (1973),  Gay  and  Preece  (1975,  1977)],  it  is  found  to  be  less  efficient. 
Finally,  in  a  recent  Naval  Postgraduate  School  (NPS)  thesis  [Ellis  (1993)],  a  pipeline 
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method  formalism  based  on  the  Cholesky  method  for  matrix  reduction  has  been 
investigated.  Although  the  mathematical  analysis  of  the  scheme  is  fairly  efficient,  the 
physical  modeling  did  not  encompass  conditions  for  other  than  laminar  flows. 

The  main  purpose  of  this  report  is  to  extend  the  range  of  flow  to  include  transition 
and  turbulent  flow  regimes.  As  this  is  primarily  a  feasibility  study,  various  facilitating 
assumptions  were  made.  First,  it  is  assumed  that  English  units  are  preferred  for 
input/output.  Second,  all  pipes  are  assumed  to  have  uniform  circular  cross-sections. 

2.      Physical  Modeling 

For  incompressible  fluids,  two  of  Bernoulli's  rules  [Lindeburg  (1982)]  for  conservation 
of  mass  and  energy  are,  the  continuity  equation 

Q  =  AlV1  =  A2V2  (1) 

and  Bernoulli's  law 

/t/2         T/2\ 

^  =  P[  2  ~    l'  +  (P2  -  Pi)  +  pg(z2  -  zi)  +  pgh'  (2) 

where  the  subscript  2/1  refers  to  positions  2/1  and  the  symbol  lists  given  in  Section 
10  defines  the  other  relevant  parameters.  The  head  loss,  as  defined  by  the  Darcy- 
Weisbach  equation,  is  given  by 

where  L  is  the  pipe  length,  d  is  the  hydraulic  diameter,  and  /  is  the  friction  factor. 
For  circular  cross-section  pipes  the  hydraulic  diameter  and  the  physical  diameter  are 
identical.  Using  similitude  analysis  [Vennard  (1961)],  it  is  possible  to  show  that  the 
friction  factor  depends  only  on  two  dimensionless  parameters,  the  Reynolds  number 

Re=^  (4) 

P 
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where  \i  is  the  dynamic  viscosity,  and  the  relative  roughness 


e 


«-5  (») 

where  e  is  the  roughness  of  the  pipe. 

A  sourceless  section  of  pipe  which  is  horizontal  [z2  =  Zi),  and  with  constant 
cross-section  so  that  V2  =  V\,  will,  from  eq  (2)  and  eq  (3),  then  have  a  pressure  drop 

Pl  -  p2  =  h  =  pgh  =  (6) 


or  after  use  of  eq  (1) 


*-'ffr>  " 


This  suggests  the  useful  electrical  analogy  with  Ohm's  law  where  pressure  drop  cor- 
responds to  voltage  and  fluid  flow  corresponds  to  current.  It  follows  that  consistent 
with  the  analogy,  the  effective  resistance  of  the  pipe  is  then  given  by 

«-4S 

where  the  units  of  this  resistance  are  defined  by  the  ratio  of  pressure  to  flow. 

It  is  observed  that  the  problem  of  predicting  fluid  flow  in  pipes  now  reduces  to 
accurate  characterization  of  the  friction  factor.  For  the  laminar  flow  domain  the 
analytically  derivable  [Vennard  (1961)]  result 

'-£ 

agrees  very  well  with  measurement  out  to  Re  <  2100  which  is  the  well  known  limit  in 
pipes,  ducts,  tubes,  and  channels  for  laminar  flow.  After  substitution  of  eq  (9)  into 
eq  (8)  it  is  found  that: 
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and  the  corresponding  relationship  for  head  loss  is 

This  is  in  agreement  with  an  expression  cited  by  Bending  and  Hutchison  (1993).  The 
subscript  I  has  been  used  to  indicate  the  laminar  flow  domain.  Note  from  eq  (10) 
that,  as  expected,  the  pipe  resistance  is  independent  of  the  flow  rate  in  the  laminar 
regime. 

In  order  to  easily  characterize  the  degree  to  which  the  fluid  flow  is  turbulent,  the 
head  loss  can  be  expressed  as 

h  =  hexrj>  (12) 

where  the  turbulence  factor,  denoted  t/>,  is  a  linear  gauge  of  the  deviation  of  the 
system  from  laminar  flow.  After  noting  that  eq  (7)  can  be  reexpressed  as 

32L/1         fpVd 

h  =  -^xQx-^r  (13) 

it  can  be  concluded  from  eq  (4)  that,  in  agreement  with  the  definition  of  the  Reynolds 
number,  the  turbulence  factors  satisfy 

*  =  /ff  (") 

As  seen  from  eq  (9)  rp  =  1  for  laminar  flow.  In  the  next  section  it  will  be  shown  that 
for  transition  and  turbulent  flow  ij>  >  1. 

The  pipe  resistance  can  now  be  expressed  as 

R  =  RexiP  (15) 

where  as  previously  noted,  Ri  does  not  depend  on  the  flow  rate. 


3.      Friction  Factor  Estimation  for  Nonlaminar 
Fluid  Flow 

For  purposes  of  the  development  here  it  is  necessary  to  introduce  two  distinct  non- 
laminar  regions.  First,  the  transition  region  satisfying 

Re0  >  Re  >  2100  (16) 

where  Re0  is  usually  cited  typically  between  3500  and  4000  [Vennard  (1961)]  and  for 
the  turbulent  region 

Re>Re0  (17) 

One  of  the  most  commonly  employed  empirical  formulas  for  the  turbulent  flow  domain 
defined  by  eq  (17)  is  attributed  to  Colebrook  (1939) 


^  =  -2.og 


e  2.51 

+ 


(18) 


3.7d      ReVJ 

Although,  as  seen  in  eq  (18),  the  dependence  of  the  friction  factor  on  Reynolds  number 
is  implicit,  the  expression  will  converge  rapidly  upon  iterative  substitution. 

To  date  the  most  viable  predictors  for  the  transition  region,  see  eq  (16),  are 
empirical.  The  Stanton  diagram  [Vennard  (1961)]  shows  smoothly  varying  curves 
which  are  based  on  measurements  in  the  transition  region.  The  Moody  diagram, 
which  only  plots  expressions  for  eq  (9),  for  Re  <  2100,  and  eq  (18),  for  Re  >  Re0,  is 
not  useful  for  the  transition  region.  Bending  and  Hutchinson  (1973)  suggested,  in  a 
computer  analysis  of  piping  networks,  that  a  linear  interpolation  between  the  laminar 
and  the  turbulent  domain  be  made  to  provide  friction  data  in  the  transition  region. 
For  the  network  analysis  scheme  to  be  described  herein,  a  linear  approximation  for 
the  transition  region  was  found  to  be  insufficiently  accurate  upon  numerical  testing. 
In  particular,  the  iterative  program  with  the  linear  approximation  would  either  not 


converge  or  converge  very  slowly.  The  problem  in  this  scheme  was  the  dramatic 
discontinuity  in  the  first  derivative  at  the  edges  of  the  transition  region.  A  greatly 
improved  modeling  scheme,  in  terms  of  convergence,  is  given  by 

'tt  (Re -2100) 


_Re_ 
*  :=  2100  + 


/(Re0)  - 


sin 


2  (Re0  -  2100) 


(19) 


2100J 

where  /(Re0)  was  calculated  via  the  Colebrook  model  eq  (18).  Note  that,  like  the 
linear  approximation  model,  continuity  is  preserved  at  the  edges  of  the  transition 
region.  Because,  as  seen  from  eq  (19),  the  slope  at  the  transition  edges  is  zero,  the 
discontinuity  in  the  slope  is  much  less  dramatic  and  the  convergence  was  found  to  be 
much  more  rapid.  Figure  1  is  a  graphical  representation  for  these  curves  in  standard 
log-log  format  for  various  relative  roughness  factors. 

It  is  worth  noting  that  the  work  by  Churchill  (1977)  agreed  very  well  with  the 
scheme  just  described  for  Re  >  2100  and  provided  an  explicit  formula  for  the  friction 
factor,  /.  However,  for  Re  <  2100  (laminar  flow)  the  expression  deviated  significantly 
from  the  classical  expression  of  eq  (9).  Hence  it  was  not  used. 

The  argument  supporting  the  claim  that  the  turbulence  factor  eq  (14)  is  generally 
greater  than  unity  is  loosely  based  on  the  following  observations.  First,  in  1913  Blasius 
proposed  an  empirical  relationship  for  smooth  pipe,  i.e.,  e  =  0  in  the  turbulent  flow 
regime  [Vennard  (1961)] 

/«nooth  pipe  =  ^        3000  <  Re  <  100, 000  (20) 

Second,  the  construction  algorithm  for  predicting  /  in  the  transition  region,  eq  (19) 

requires  that 

64 

/transition  _  /laminar  —   p  \~^J 

This  leads  to  the  inequality  that  is  established  from  a  comparison  of  eq  (20),  eq  (21), 


and  Fig  1 
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/rough  pipe  _  /smooth  pipe  _  J  laminar  —    p  \^^) 


and  is  supported  by  analysis  and  measurement  [Vennard  (1961)].  The  inequality, 
rp  >  1,  then  follows  directly  from  the  definition  of  the  turbulence  factor,  xp,  given  by 
eq  (14). 

The  turbulence  factor  has  been  introduced  into  this  report  for  two  reasons.  Al- 
though the  friction  factor  is  a  linear  indicator  for  the  system  response  of  the  fluid 
dynamics,  it  is  not  a  linear  indicator  for  the  deviation  of  the  system  from  laminar 
behavior.  Note  that  the  quantity  (Re  —  2100)  is  a  nonlinear  indicator  for  the  devia- 
tion from  laminar  behavior.  The  purpose  of  this  work  has  been  to  extend  the  Ellis 
(1993)  effort  to  encompass  nonlaminar  regimes.  The  turbulence  factor  will  readily 
demonstrate  the  need  for  this  extension.  Second,  it  turns  out  that  the  laminar  pipe 
resistance,  /?*,  can  be  computed  at  the  onset  of  the  problem  and  only  the  turbulence 
factors  need  to  be  updated  in  the  iterative  calculations  of  eq  (15)  rather  than  to  com- 
pute eq  (8).  In  brief,  the  turbulence  factors  provide  some  computational  advantages 
which,  for  large  piping  networks,  would  be  significant. 

4.      Basic  Circuit  Modeling 

Following  the  electrical  analogy  for  an  arbitrary  kth  pipe,  a  pressure  difference,  Ap*, 
across  a  pressure  source,  Apak,  in  series  with  a  pipe  resistance,  r*,  with  a  flow  </*, 
satisfies 

Ap*  =  qkrk  -  Apak  (23) 

consistent  with  the  convention  established  by  Fig  2.  The  corresponding  relation  in 
terms  of  a  current  source  is  given  by 

qk  =  q„k  +  Vk&Pk  (24) 
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Figure  1:     Friction  factor  model. 
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Figure  2:     Pressure  source  circuit  model. 


where  y^  =  \/rk.  Figures  2  and  3  define  the  nomenclature  and  the  topology.  Equation 
(24)  can  be  compactly  represented  for  all  branches  or  pipes  in  matrix  format 


Q  =  Q3+YAP 


(25) 


where  it  is  seen  that  eq  (24)  and  eq  (25)  have  unresolved  unknowns  associated  with 
flow  rate  and  pressure.  The  condition  of  continuity  at  the  nodes 


J2^j  =  Q        ;  =  1,2,  ...nT 
k 


(26) 


where  n  j  is  the  total  number  of  nodes  in  the  network,  provides  enough  information  to 
resolve  the  flow  rates  and  pressure  drops.  A  formulation  for  this  process  is  presented 
in  the  next  section. 
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Figure  3:     Flow  source  circuit  model. 

5.     Mathematical  Formulation 

The  mathematical  formalism  is  presented  via  an  illustration  taken  from  the  Ellis 
(1993)  thesis.  The  piping  network  shown  in  Fig  4a  has  five  nodes  and  six  connecting 
branches.  It  is  also  seen  that  the  branch  between  nodes  1  and  5  contains  a  pump. 
This  pump  provides  a  source  for  fluid  flow  at  its  high  pressure  side  shown  at  node-1 
with  an  arrow.  The  low  pressure  side  is  the  corresponding  sink.  The  corresponding 
flow  graph  is  shown  in  Fig  4b.  Each  branch  now  has  an  arbitrarily  chosen  orientation 
arrow  which  defines  the  convention  for  positive  fluid  flow.  Note  that  consistent  with 
the  flow  direction  shown  for  branch  1,  the  source  shown  on  Fig  4a  would  be  assigned 
a  negative  numerical  vale.  The  first  step  in  the  mathematical  formalism  is  to  cre- 
ate both  a  network  representation  for  the  fluid  pipeline  network  and  a  corresponding 
flow  graph.  The  flow  rates  in  the  b  branches,  in  this  case  6,  can  be  represented  as  a  6x  1 
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Figure  4:     a)  Six  branch  pipeline  network;  b)  Associated  flow  graph  for  a)  [15]. 
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column  vector 


Q  = 


9i 

?2 
93 
94 

9fl  J 


(27) 


The  elements  of  the  augmented  node-branch  incidence  matrix,  c° '•,  can  be  defined 


as 


c*     =        1         Branch  j  leaves  node  i 


c"j    =    —  1         Branch  j  enters  node  i 
=       0         otherwise 


'«j 


(28a) 
(28b) 
(28c) 


where  i  =  1, 2, . . . ,  nj  nodes  and  j  =  1, 2, . . . ,  6  branches.  This  leads  to: 


Ca  = 


110  0  0  0 
0-11100 
0  0-1010 
0  0  0-1-1  1 
-10       0       0       0-1 


(29) 


where,  in  general,  Ca  has  6  columns  and  nj  rows.  It  is  apparent  that  for  each  column, 
which  is  associated  with  a  node, 


Ec^  =  0 


(30) 


or  more  formally 


CaQ  =  0  (31) 

consistent  with  flow  rate  conservation  at  all  nodes  (see  eq  (26)).  In  order  to  reduce 
the  order  of  the  problem  by  1  it  is  useful  to  define  one  of  the  nodes  as  the  ground  or 
datum  reference.  In  principle  it  should  not  matter  which  node  is  taken  as  this  datum 
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node,  however,  in  the  coded  implementation  it  is  assumed  that  the  circuit  nodes 
are  numbered  such  that  the  datum  node  is  the  node  bearing  the  highest  numerical 
designator.  For  example,  as  seen  on  Figs  4a,  b,  the  node-5  is  the  datum  node.  Because 
the  assignment  of  nodes  is  independent  of  the  assignment  of  branches  and  is  arbitrary, 
there  is  no  sacrifice  in  the  generality  of  the  method.  Following  this  point,  the  node- 
branch  incidence  matrix,  C,  can  be  defined  according  to  the  rules  dictated  by  eq  (27) 
with  the  index  i  satisfies  the  condition  i  =  1,2, ...  ,n  where 


n  =  nj  —  1 


Thus,  for  this  example,  the  node-branch  incidence  matrix  is 


(32) 


C  = 


110  0  0  0 
0-1  1  1  0  0 
0  0-1010 
0       0       0-1-11 


(33) 


In  general  C  has  6  columns  and  n  rows.  It  follows  from  eq  (30)  that 


CQ  =  0 


(34) 


which  forces  a  condition  of  flow  rate  conservation  or  continuity  at  all  nodes. 

The  node  pressures,  defined  with  respect  to  the  reference  datum  node,  can  be 
represented  by  the  column  vector 


P  = 


Pi 

V2 
Pn   J 


(35) 


which,  for  the  example  under  consideration,  the  node  pressure  vector  P  is  a  4  x  1 
column  vector  since  n  =  4.  The  node  branch  incidence  matrix  can  be  used  to  predict 
the  head  losses  from  the  rule 

H  =  CTP  (36) 
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where  the  branch  head  loss  vector,  represented  as: 


H  = 


hi 
h2 


L  riB 


(37) 


can  be  calculated  using  eqs  (11)  and  (12).  For  the  example  being  examined,  H  is  a 
6x1  column  vector,  corresponding  to  the  number  of  branches. 

In  Section  4  the  pressure  drop  vector  was  introduced  without  any  assumptions 
regarding  the  physical  modeling  introduced  in  Section  2.  Under  the  practical  as- 
sumptions made  leading  to  eq  (6),  it  is  correct  to  let 


AP-^H 


in  eq  (25),  which  leads  to 


Q  =  Q5  +  YH 


(38) 


(39) 


Premultiplication  of  this  equation  by  C  leads  to: 


0  =  CQs  +  CYH 


(40) 


where  the  left  hand  side  of  eq  (40)  is  predicted  from  eq  (34).  After  substitution  of  eq 
(36)  into  eq  (40)  and  solving  for  P,  it  follows  that 


-i< 


P  =  Y-Q 


(41) 


where 


and 


Q  =  -CQ, 


Yn  =  CYCJ 


(42a) 


(42b) 
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Given  the  matrix  Yn  and  the  source  vector  Q5,  any  number  of  matrix  inverting 
schemes  [Hamming  (1962)]  could  be  employed  to  evaluate  eq  (40).  As  suggested  in  the 
Ellis  (1993)  thesis  the  very  efficient  Cholesky  reduction  algorithm  is  applicable  here 
because  the  matrix  Yn  is,  not  only  symmetric,  but  positive  definite.  The  Cholesky 
reduction  into  upper  and  lower  triangular  matrices  permits  a  rapid  matrix  inversion 
using  Gaussian  elimination  as  described  in  Appendix  A. 

Once  the  node  pressures  given  by  eq  (40)  are  obtained  the  branch  head  losses  and 
flow  rates  are  easily  calculated  from  eq  (36)  and  eq  (39),  respectively.  Note  that,  once 
the  flow  rate  vector  is  obtained,  the  flow  rate  velocities  can  be  computed  according 
to  the  rule 


V  = 


vb 


=  -Q 


(43) 


in  agreement  with  eq  (1).  The  corresponding  velocity  magnitudes,  (|i>i|,  |t>2|,  •  •  • ,  \vb\) 
can  then  be  applied  with  eq  (4)  to  calculate  the  Reynolds  numbers,  the  friction  factors, 
and  the  turbulence  factors  using  eq  (14).  The  pipe  resistances  are  calculated  from  eq 
(15)  and  the  corresponding  admittances,  2/1,1/2,-  •  •  ,2/6,  follow  from 


Y  = 


0 

0 
l/r2 

0     •• 

•      0 

0 

0 

'•. 

0 

(44) 


:  ■•■      0 

0         0       •••      0      l/r6 

If  the  current  turbulence  factors  differ  significantly  from  the  previous  turbulence  fac- 
tors according  to  the  rule 

MAXERR  >  max  ^(previous)  -  V>it ( current )|    :      {branches  k  =  1,6}         (45) 

the  process  is  repeated.  The  Fortran  input  parameter  MAXERR  defines  the  maximum 
deviation  in  the  turbulence  factors  between  iterations.  It  is  assumed  initially,  in  the 
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first  step  of  the  iterative  scheme,  that  the  flow  is  in  the  laminar  regime,  that  is,  all 
turbulence  factors  =  1. 

6.      Input/Output:  Parameters,  Conventions, 
Procedures 

It  follows  for  the  definition  of  specific  weight  that 

7  =  99  (46) 

where  g,  the  gravitational  constant,  is  32  ft/sec2.  Using  eq  (46)  the  computer  program 
computes  the  density  which  is  needed  for  calculation  of  the  Reynolds  numbers  eq  (8). 
Figure  5  defines  the  convention  for  distinguishing  the  entry  (F)  and  exit  (T)  nodes 
of  a  particular  branch.  The  arrow  defines  the  direction  of  positive  flow. 


F  T 


Figure  5:     Convention  defining  nodes. 

The  required  inputs  to  the  FORTRAN  program  code  flow. for  provided  in  Ap- 
pendix A  are  displayed  in  Table  1. 

The  calculations  are  done  in  formal  English  units  of  (ft,  Lb,  sec).  However,  the 
inputs  and  outputs  are,  when  appropriate,  given  in  conventional  units.  Conversion 
factors  for  these  cases  are  shown  in  Table  2. 
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Table  1:     Input  Variables 


Symbol 

Section 

Units 

rieo 

§3 



typical  3500 

MAXERR 

§5 

— 

typical  0.05 

HLCF 

§2 

— 

typical  60.0 

TIJ 

§5 

— 

depends  on  network 

b 

§5 

— 

depends  on  network 

V 

§2 

lO"5  lb-sec/ft3 

typical  H20,  3.75 

w 

§6 

lb/ft3 

typical  H20,  64.0 

node  specification  —  .F,  —  T 

§6 

— 

depends  on  pipe 

L 

§2 

ft 

depends  on  pipe 

d 

§2 

inches 

depends  on  pipe 

e 

§2 

mil 

depends  on  pipe 

source  specification  Qa 

§5 

gal/min 

depends  on  pump(s) 

Table  2:     Conversion  Factors 


Formal 

Conventional 

Conversion 

Description 

Symbol 

Units 

Units 

Factor 

Flow  rate 

Q 

ft3/sec 

gal/min 

450 

Pressure 

p 

lb/ft2 

psi 

1/144 

Roughness 

e 

ft 

mil 

12,000 

In  order  to  facilitate  the  process  of  batch  input,  the  code  has  been  designed  to 
read  the  input  file  flowinp.dat  rather  than  require  the  user  to  laboriously  submit 
the  data  in  the  interactive  mode.  In  this  way  a  maximum  amount  of  flexibility  can 
be  incorporated  into  the  program  without  creating  a  tedious  data  entry  procedure  for 
the  user  of  the  code.  Thus,  the  "user  friendly'1  goal  has  been  kept  in  mind. 
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The  outputs  from  the  code  flow. for,  which  are  listed  by  branch  number  (BR), 
are  displayed  in  Table  3. 

Where  appropriate  both  text  symbol  and  fortran  symbol  have  been  specified  in 
Table  4.  For  example,  PF  and  PT  refer  to  node  pressures  at  entry  and  exit,  respec- 
tively. In  some  practical  cases  it  may  be  desirable  to  employ  the  equivalent  of  an 
"ideal"  constant  flow  rate  pump,  that  is,  qk  =  qak  in  eq  (24).  This  is  readily  handled 
by  forcing  yk  — *  0  or  equivalently  rk  — *  oo.  Because  of  the  head  loss  coupling  factor 
rule-of-thumb  introduced  at  the  end  of  Section  3,  the  most  convenient  way  to  force 
the  pump  source  to  behave  ideally  is  to  let  dk  —*  0.  In  practice,  because  of  numerical 
instabilities,  if  yk  =  0,  it  is  recommended  that  the  user  set 


dk  =  10      x  min{dt}         i  =  1,2, ...  ,6        but  i  ^  k 


(47) 


A  hardcopy  of  the  input/output  for  the  piping  networks  tested  is  provided  in 
Appendix  B.  A  more  detailed  description  of  these  cases  is  covered  in  the  next  section. 

Table  3:     Output  Variable 


Symbol 

Section 

Units 

BR 

§5 

_ 

Q 

§5 

gal/min 

PF 

§6 

10~3  psi     (mpsi) 

PT 

§6 

10~3  psi     (mpsi) 

HEAD  =  H 

§5 

10~3  psi     (mpsi) 

ERR 

see  RHS  eq  (44) 

— 

Re 

§2 

— 

TBF  =  V 

§2 

— 

Y 

§5 

(gal/min) /mpsi 
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7.      Examples 

General  Observations 

Several  general  observations  can  be  easily  confirmed  from  the  data  (see  Appendix  B 
for  all  the  case  examples  discussed  in  this  section).  Specifically, 

PF  -  PT  =  HEAD         (pF  -pT  =  hk) 

and 

Q  =  Qa  +  Y  HEAD        (qk  =  qak  +  ykhk) 

In  all  the  examples  the  sink  node  of  the  pump  was  arbitrarily  taken  as  the  datum 
node  for  convenience  of  interpretation.  The  node  pressures  under  columns  PF  and 
PT  are  as  expected,  strictly  positive.  In  all  three  examples  some  or  all  of  the  pipeline 
segments  were  operating  in  the  nonlaminar  regime.  This  is  evident  by  checking  the 
TBF  columns  of  each  example  in  Appendix  B.  The  associated  Reynolds  number  is 
also  provided.  The  branch  errors  were  calculated  according  to  eq  (45).  It  is  easily 
checked  that  the  MAXERR  parameter,  line  2  of  the  input  set,  put  a  ceiling  on  these 
branch  error  values. 

CASE  A  —  A  four  branch  pipeline  network  with  an  ideal  source  (see  Fig  6a  and  6b). 
This  case  illustrates  the  generation  of  an  ideal  flow  source.  Note,  as  discussed,  this 
condition  can  be  obtained  by  making  the  pipe  diameter  of  the  source  branch  relatively 
small.  The  associated  Reynolds  number  and  turbulence  factors  will  be  artificially  high 
and  may  not,  as  represented  in  the  line  1  output  for  Case  A,  be  within  the  defined 
format. 
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WW\A 


Figure  6:     a)  4  branch  with  ideal  source;  b)  4  branch  flow  graph  for  a). 

CASE  B  —  A  six  branch  pipeline  network  with  nonideal  source  (see  Fig  4a  and  4b). 
Note,  as  seen  here,  the  effect  of  the  nonzero.admittance  of  node  1  is  to  impede  the  flow 
of  the  source  into  the  rest  of  the  pipeline  circuit.  Also,  the  source  is  required  to  be 
negative  because  the  pump  is  physically  driving  the  flow  opposite  to  the  convention 
defined  by  the  arrow  of  branch  #1. 
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CASE  C  —  A  12  branch  piping  network  for  shoebox  cooling  rack  (see  Fig  7).  Note 
that  several  of  the  branches  have  negative  Q  values.  The  physical  interpretation  is 
that  these  cases  have  a  fluid  flow  opposite  to  the  convention  defined  on  Figure  7. 


§  W  4 

Figure  7:     12  branch  shoebox  cooling  rack. 


8.      Future  Enhancements  in  the  Modeling 

Although  the  work  described  herein  extends  the  Ellis  (1993)  thesis  by  including  the 
nonlaminar  flow  regime,  it  still  should  be  considered  as  preliminary.  For  example: 

•  The  network  formalism  discussed  herein  should  be  compared '  with  the  bench- 
mark Hardy  Cross  method  [Lindeburg  {1987)]  as  well  as  other  more  recently 
proposed  schemes  [Carnahan  and  Christensen  (1972),  Bending  and  Hutchinson 
(1973),  Gay  and  Preece  (1975,  1977)]  in  order  to  assess  the  relative  computa- 
tional efficiencies. 

•  The  next  major  modeling  development  in  the  algorithm  would  be  the  addition 
of  the  capability  of  evaluating  the  node  to  node  heat  transfers. 
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•  Using  the  theory  for  compressible  fluid  flow  [Vennard  (1961)]  required,  for  ex- 
ample, when  the  cooling  fluid  is  a  gas,  the  developed  scheme  could  be  extended 
to  allow  for  compressible  fluids.  Due  to  the  sensitivity  of  gas  properties  to  tem- 
perature this  development  would  have  to  include  effects  due  to  heat  transfer. 

•  An  alternate  but  similar  mathematical  formalism  could  be  derived  for  pumps 
best  characterized  as  constant  pressure  sources.  It  is  certainly  possible  to 
"Theveninize"  the  source  in  the  context  of  the  present  mathematical  formal- 
ism. However,  it  should  be  noted  that  because  of  the  nonlinearity  inherent  in 
the  model,  this  Theveninization  would  need  to  be  repeated  iteratively,  decreas- 
ing the  efficiency  of  the  process. 

•  As  previously  noted  in  Section  3,  a  significant  improvement  in  the  rate  of  con- 
vergence of  the  method  was  obtained  by  reducing  the  discontinuity  in  the  fric- 
tion factor  derivatives  at  the  transition  edges.  It  should  be  possible  to  apply 
mathematical  adjustments  on  eq  (19),  for  example,  using  a  technique  known  as 
Hermiteor  "osculating"  interpolation  [Hamming  (1962)],  in  order  to  completely 
eliminate  the  slope  discontinuity.  The  potential  gain  in  the  rate  of  convergence 
could  be  significant. 

•  Various  mundane,  but  no  doubt  practical  embellishments  could  also  be  consid- 
ered. For  example,  an  option  of  metric  system  units  could  be  added.  Also,  an 
option  to  specify  pipe  bends,  vertical  elevation,  and  cross-section  types  could  be 
added.  If  commercial  viability  is  of  concern,  the  program  could  be  dovetailed 
with  CAD  software. 

•  The  required  background  investigation  revealed  that  despite  a  long  history  of 
developments  on  modeling  fluid  flow  in  pipes,  only  recently  have  mathematical 
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methods  involving  the  theory  of  chaos  begun  to  unravel  and  predict  the  prop- 
erties of  nonlaminar  fluid  flow.  A  small-scale  physical  model  which  predicts  all 
the  measured  properties  on  the  Stanton  diagram  is  at  this  date  unavailable. 

9.      Conclusion 

A  mathematical  model  for  predicting  laminar  and  nonlaminar  flows  in  pipeline  net- 
works has  been  proposed  and  successfully  tested.  This  scheme,  upon  iterative  appli- 
cation, has  been  found  to  converge  to  an  approximate  solution.  The  computational 
efficiency  of  this  convergence  process  was  found  to  be  highly  dependent  on  the  fric- 
tion factor  modeling  in  the  transition  region  between  laminar  and  turbulent  behavior. 
In  particular,  a  significant  improvement  in  the  convergence  rate  was  obtained  by  re- 
placing the  previously  proposed  linear  model  with  a  nonlinear  model  having  a  less 
dramatic  first  derivative  discontinuity.  A  linear  gauge  for  the  deviation  from  laminar 
behavior,  referred  to  as  the  turbulence  factor,  was  introduced  in  this  report  in  order 
to  facilitate  the  computational  task. 
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10:      Nomenclature 


Roman  Letter  Symbols 


Units 


A 

pipe  cross-sectional  area 

ft2 

d 

pipe  diameter 

ft 

e 

pipe  roughness 

ft 

f 

friction  factor 

— 

9 

gravitational  acceleration 

ft /sec2 

h 

head  loss  due  to  friction 

lb/ft2 

P 

pressure 

lb/ft2 

Q 

flow  rate 

ft3/sec 

R 

pipe  resistance 

Lb-sec/ft5 

Re 

Reynolds  number 

— 

V 

velocity 

ft /sec 

V>3 

time- rate- change  of  work  output  (pump) 

lb-ft/sec 

Y 

pipe  admittance 

ft3/(lb-sec) 

z 

elevation 

ft 

Greek  Letter  Symbols 


7  specific  weight 

t  pipe  roughness  ratio 

fi  viscosity 

p  density 


lb/ft3 

lb-sec/ft3 
slug/ft3 
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Appendix  A:      Cholesky  Reduction 

This  appendix  provides  a  more  complete  description  on  the  Cholesky  reduction  algo- 
rithm [Kraus  (1992)].  Positive  definite,  symmetric  matrices  can  be  factored  according 
to  the  rule 

M  =  LLT  (Al) 

where  L  is  an  upper  triangular  matrix.  Because  of  the  applicability  of  this  theorem 
[Ellis  (1993),  Kraus  (1987)]  to  the  network  matrix  Yn  eq  (41b)  it  follows  that 

Yn  =  LLT  (A2) 

and  for  eq  (40) 

LLTP  =  Q  (A3) 

The  advantage  of  the  factorization  eq  (Al)  is  now  apparent  since  eq  (A3)  can  be 
broken  into  two  parts:  First, 

LJ  =  Q  (A4a) 

and  second, 

LTP  =  J  (A4b) 

which  are  easily  solved  in  succession  via  Gaussian  elimination  to  determine  first  the 
column  vector  J  as  an  intermediate  step  then  the  derived  node  pressures  P. 
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Appendix  B:      Input  /Output  Hardcopy 


****************   INPUT  FORMAT  ************************* 

REYNOLDS'  NUMBER  FOR  RIGHT  EDGE  OF  TRANSITION  REGION 

ERROR  TOLERANCE  IN  CALCULATION 

THE  HEADLOSS  COUPLING  FACTOR  (TYPICAL  NUMBER  60  ) 

THE  NUMBER  OF  NODES 

THE  NUMBER  OF  BRANCHES 

THE  VISCOSITY  X  10  -5  LB-SEC/FT~3 

THE  DENSITY  LB/FT  ~3 

STARTING  NODE,  ENDING  NODE,  LENGTH(FT) .DIAMETER(IN) ,  ROUGHNESS (MILS) 

NOTE  FOR  IDEAL  CURRENT  SOURCE  SET  DIAM=.001X  SMALLEST 
THE  NUMBER  OF  SOURCES 
SOURCE  STRENGTH  GAL/MIN 
SOURCE  BRANCH 

************************************************************ 
INPUT  (CASE  A  4  BRANCHES  IDEAL  SOURCE) 
3500.0 
0.005 
60.0 
4 
4 

3.75 
64. 

4.1,  1.0, .001,1.0 

1.2,  1.0,1.0,2.0 


2,3,  1 

.0,1. 0,( 

).5 

3.4,  1 

1 

6.0 

1 

.0,1.0,0.6 

OUTPUT 

CASE  A 

BR 

Q 

PF 

PT 

HEAD 

ERR 

Re 

TBF 

Y 

1 

6.00 

.0 

282.1 

-282.1 

.0039 

******* 

******* 

.0000 

2 

6.00 

282.1 

183.8 

98.2 

.0032 

10806.2 

5.6 

.0609 

3 

6.00 

183.8 

92.1 

91.7 

.0032 

10806.2 

5.2 

.0652 

4 

6.00 

92.1 

.0 

92.1 

.0032 

10806.2 

5.3 

.0649 

*********************************************************************** 

INPUT  (CASE  B  6  BRANCH  PIPELINE  NETWORK) 

3500.0 
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0.01 

60.0 

5 

6 

3.75 

64. 

1,5,  1.0,1.0,1.0 

1.2,  1.0,1.0,2.0 

2.3,  1.0,1.0,0.5 

2.4,  1.0,1.0,0.6 

3.4,  1.0,1.0,0.6 

4.5,  1.0,1.0,5.0 
1 

-6.0 
1 


OUTPUT  CASE  B 


BR 


Q 
-1.75 
1.75 

.58 
1.17 

.58 
1.75 


PF 

25.1 

25.1 

14.6 

14.6 

12.9 

11.2 


PT 

.0 
14.6 
12.9 
11.2 
11.2 

.0 


HEAD 


25 
10 

1 

3 

1 

11, 


ERR 
.0073 
.0075 
.0000 
.0000 
.0000 
.0078 


Re 
3149.7 
3149.7 
1049.9 
2099.8 
1049.9 
3149.7 


TBF 
2.0 
2.1 
1.0 
1.0 
1.0 
2.2 


Y 

.1679 
.1645 
.3409 
.3409 
.3409 
.1554 


*********************** **** **********************  ******************** 

CASE  C  (  12  BRANCH  SHOEBOX  COOLING  RACK) 

3500.0 

0.01 

60.0 

8 

12 

3.75 

64. 

8,2, 

2,3, 


0,1.0,1.0 
0,1.0,2.0 
3,4,  1.0,1.0,0.5 
4,8,  1.0,1.0,0.6 
5,8,  1.0,1.0,0.6 

2.6,  1.0,1.0,5.0 

7.3,  1.0,1.0,1.0 

1.4,  1.0,1.0,2.0 

6.5,  1.0,1.0,0.5 

6.7,  1.0,1.0,0.6 
7,1,  1.0,1.0,0.6 
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5,1.  1 
1 

6.0 
1 

.0,1.0,5 

.0 

OUTPUT 

CASE  C 

BR 

Q 

PF 

PT 

HEAD 

ERR 

Re 

TBF 

Y 

1 

3.48 

.0 

24.6 

-24.6 

.0341 

6273.0 

3.4 

.0990 

2 

1.76 

24.6 

14.5 

10.1 

.0393 

3168.5 

2.0 

.1681 

3 

1.33 

14.5 

9.6 

4.9 

.0084 

2394.5 

1.3 

.2685 

4 

1.74 

9.6 

.0 

9.6 

.0371 

3138.2 

1.9 

.1750 

5 

1.74 

9.6 

.0 

9.6 

.0366 

3134.8 

1.9 

.1752 

6 

1.72 

24.6 

14.4 

10.2 

.0386 

3104.5 

2.1 

.1635 

7 

-.43 

13.2 

14.5 

-1.3 

.0000 

773.9 

1.0 

.3409 

8 

.41 

10.8 

9.6 

1.2 

.0000 

743.6 

1.0 

.3409 

9 

1.32 

14.4 

9.6 

4.8 

.0072 

2381.3 

1.3 

.2710 

10 

.40 

14.4 

13.2 

1.2 

.0000 

723.2 

1.0 

.3409 

11 

.83 

13.2 

10.8 

2.4 

.0000 

1497.1 

1.0 

.3409 

12 

-.42 

9.6 

10.8 

-1.2 

.0000 

753.5 

1.0 

.3409 

28 


Appendix  C:      Computer  Code 


$debug 

$LIST 

c  program  for  CALCULATING  PRESSURES  AND  FLOW  RATES  IN  PIPILINE  NETWORKS 

C  ******************************************************************* 

C  PROGRAM  flow. FOR  WRITTEN  BY  PROF.  RON  J  PIEPER 

C  NPS  408  6562101 

C  AUG.  30,  1994 

C  SEE  FILE  flowINP.DAT  FOR  INPUT 

C  SEE  FILE  flowOUT.DAT  FOR  OUTPUT 

C  SEE  FILE  flowDIA.DAT  FOR  DIAGNOSTICS  ON  INPUT 

C  DESIGNED  TO  HANDLE  UP  TO  40  BRANCHES 

C  CAN  BE  ADJUSTED  BY  INCREASING  DIMENSION 

C  OF  THE  ARRAYS 

C  ******************************************************************** 

c  mu  fluid  viscosity  (lb-sec/ft"2) 

c  rho  fluid  density   (lb/ft~3) 

c  g  acceleration  due  to  gravity  (ft/sec~2) 

c  nu  specific  gravity  rhoxg 

c  pi  pi 

c  ren#  reynold's  •  (rho*vel*d/mu) 

c  f  friction  factor 

c  vel  velocity  of  fluid 

c  pres  pressure  of  the  fluid 

c  *****************pipe  variables  ************************************** 

c  ell  pipe  length(s) 

c  d  pipe  diameter(s) 

c  e  roughness 

c  rat  roughness  ratio 

C  HCLF  HEADLOSS  COUPLING  FACTOR 

c  **********  p  FACTOR   ******************************** 

C  TURB  TURBULENCE  FACTOR  =  F*2100/64  (=1.0  LAMINAR  REG) 

C  MAXERR  MAXERR  IN  TURBULENCE  FACTOR 

c  FOR  REFERENCE  ALSO  SEE  ITERATIVE  COOLBROOKE-WHITE 

C  AND  ALSO  SEE  DIRECT  CHURCHILL  eQ  (  IN  PLOT  PROGRAM) 

c  applies  for  ren#  >  REYEDGE 

C  REDGE  THE  LEFT  EDGE  OF  TRANSITION  REGION 

C  REYNOLDS#  WHERE  COOLBROOK  BEGINS 

c  ************************************************************** 

c  ******  network  constants  ************************************* 

c  N  the  numver  of  junctions  in  the  system 
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c  Nl  the  #  junction  -  datum  node  (n-1) 

c  B  the  #  branches  in  the  system 

C  S  the  #  of  sources 

c  *************************************** ********************* 

c  ASSUMPTIONS 

C  1.    C00LBR00KE-WHITE  FORMULAE  FOR  REY#>  redge 

C  2.    WORK  IN  ENGLISH  UNITS 

C  3.    FLUID  IS  INCOMPRESSIBLE 

c  4.    ALL  PIPES  ARE  CIRCULAR  IN  SHAPE 

c  5.    THE  CORRECT  SOLUTION  IS  OBTAINED  BY  ITERATION 

C  ************************************************************* 

integer  N,N1,B,QB,S,DN 

integer  NT(40),  NF(40),X,  itest 

real  mu,  rho,  g,  pi,  C(40,40),  d(40),  ell(40),  r(40) ,QS1,QS(40,1) 

REAL  RLAM(40),TURB(40),Ct(40,40),Cy(40,40),YLAM(40),  Y(40,40) 

REAL  IS(40,1),YN(40,40),P(40,1),  YH(40,1) ,V(40, 1) ,  REY(40,1) 

REAL  AREA(40) ,ISI(40,1) ,YNI(40,40) ,Q(40, 1) ,ERR(40) ,ep(40) , rat (40) 

REAL  TURBF,MAXERR,QG(40),FCROS(40) 

REAL  REDGE,  EPS, H (40, 1) ,HCLF 

CHARACTER*20  FORM 

CHARACTER*3  BR 

CHARACTER*7  NAME(ll) 

BR='BR  ' 


NAME(7)= 
NAME (4)= 
NAME(2)= 
NAME(3)= 
NAME(6)= 
NAME(1)= 
NAME (8)= 
NAME(5)= 


TBF 
HEAD' 
PF  ' 
PT  ' 
Re   ' 

Q' 
Y   » 
ERR' 


4     FORMAT  (A12) 

OPEN  (UNIT=7,FILE=' flowDIA.DAT'  .STATUS* 'NEW  ,ACCESS=' SEQUENTIAL') 
OPEN  (UNIT=8 , FILE= ' f lowINP . DAT ' ,  STATUS= ' OLD ' , ACCESS= » SEQUENTIAL ' ) 
OPEN  (UNIT=9 , FILE= ' f lowOUT . DAT ' ,  STATUS= ' NEW , ACCESS= ' SEQUENTIAL ' ) 

**************************************************************** 

g=32 . 174 
pi=3. 1415926 
READ(8,8) 
8     FORMAT  (/,/,/./././././././,/) 

C     write(*,13) 
write(7,13) 
13    format  ('INPUT  THE  REYNOLDS  •  FOR  RIGHT  EDGE  TRANSITION  ') 

30 


read(8,30)  redge 

WRITE(7,30)  REDGE 
C     write(*,23) 

write(7,23) 
23    format  ('  enter  the  maximum  error  IN  TURBULENCE  FACTOR  desired'  ) 

read (8, 30)  MAXERR 

WRITE (7, 30)  MAXERR 
C     WRITE(*,21) 
21   FORMAT  ('ENTER  THE  AVERAGE  HEADLOSS  COUPLING  FACTOR, TYPICAL  #60') 

READ(8,30)  HLCF 

write(7,21) 

write (7, 30)  HLCF 
C     write(*,10) 

write(7,10) 
10    format  ('  enter  the  number  of  nodes  in  the  rack  system  N>100  '  ) 

read(8,15)  N 

N1=N-1 
15    format (15) 

write(7,15)  N 
C     write(*,15)  N 

C     write(*,20) 

write(7,20) 
20    format  ('  enter  the  number  of  branches  in  the  rack  system  B>N  '  ) 

read(8,15)  B 
C     write(*,15)  B 

write(7,15)  B 

^^*^*****^m********** *****************  ********************************* 

C     write(*,25) 

write(7,25) 
25    format  ('  input  the  viscosity  __  X  10~-5  lb-SEC/ft"2  *  ) 

read(8,30)  mu 

WRITE (7, 30)  mu 
C     WRITE (*, 30)  mu 

mu=mu*. 00001 
30    format  (fl0.4) 
C     write(*,35) 

write(7,35) 
35    format (  '  input  the  specific  weight  lb/ft"3'  ) 

read(8,30)  sw 

rho=sw/g 
C     write(*,30)  sw 

write(7,30)  sw 
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c  Node-branch  matrix 

c 

c    matrix  initialization 

c 

do  40  i=l,Nl 
do  40  j=l,B 
40    C(i,j)=0 

do  42  i=l,B 

NF(i)=0 
42    NT(i)=0 
C 
C 

C     receive  dat  from  the  keybord  to  develop  a,d  and  ell  matrices 
c     for  circular  pasages 

c     start  loop  on  branches  to  input  node  to  node  information 
do  50  i=l,B 
write(7,55)  i 
55    format  ('  at  branch  number  ',  2x,  i5) 
C     write(*,60) 
WRITE (7, 60) 
60    formatCENTER  nf (i) ,nt(i) .length(ft) .diameter(in) ,ROUGHNESS-mil') 
read(8,*)  NF(i) ,NT(i) ,ell(i) ,d(i) ,ep(i) 
ep(i)=ep(i)*.001 
rat(i)=ep(i)/d(i) 
c   convert  from  mils  to  inches 

WRITE(7,*)  NF(i) ,NT(i) ,ell(i) ,d(i) ,ep(i) ,rat(i) 
C  WRITE(*,*)  NF(i),NT(i),ell(i),d(i),ep(i),rat(i) 
C  70    formatC  i5,i5,f8.3,f8.3,f9.5) 

c  *************  calculate  the  coolbrook  crossing  value  at  R=  3500 
c        colebrook-white  formula  page  198  of  fox,  let  f  goto  4f 
sqf=(64. 0/2100.0)**. 5 
3100        save=sqf 

sqf=1.0/(-2.0*loglO(2.51/(REDGE*save)  +  rat(i)/3.7)) 
IF  (SAVE  .EQ.  SQF)  GOTO  4000 
error=abs(save-sqf )/sqf 
if  (error  .gt.  .001)  goto  3100 
4000        fcros(I)=sqf**2 
C     WRITE(*,*)  'I,FCR0S(I),RAT(I)\  I,  FCROS(I),  RAT(I) 
50    continue 
C  SECTION  ON  SOURCE  INPUT 
c        INPUT  SOURCES,  S  OF  THEM 
C        QS  MATRIX  OF  FLOW  SOURCES 
C        QB  SOURCE  BRANCH 
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C        QS1  SOURCE  STRENGTH 
C        WRITE(*,90) 

WRITEC7.90) 
90    FORMAT (  '  INPUT  THE  NUMBER  OF  SOURCES  ') 

READ(8,92)  S 

write(7,92)  S 
C       write(*,92)  S 

92   F0RMAT(I5) 
DO  98  1=1, S 
C     VRITE(*,95)  I 
WRITE(7,95)  I 
95   FORMATC  'INPUT  THE  SOURCE  STRENGTH  FOR  THE',  15, 'TH  SOURCE') 
READ(8,30)  QS1 
write(7,30)  QS1 
C     write(*,30)  QS1 

QSl=QSl/450 
C    CONVERT  FROM  GAL/MIN  TO  FT~3/SEC  CONVERSION  FACTOR  1/450 
C    X  GAL/MIN=  (450) *-l  FT~3/SEC 

•  WRITE(7,93) 
C       WRITE(*,93) 

93    FORMATC  INPUT  THE  BRANCH  NUMBER  OF  THE  SOURCE  ') 
READ(8,15)  QB 
write(7,15)  QB 
C     writeO,15)  QB 
98    CONTINUE 

DO  100  K=1,B 
IF  (K  .EQ.  QB)  THEN 
QS(K,1)=  QS1 
ELSE 
QS(K,1)=0 
ENDIF 
100   CONTINUE 

c  Start  setting  up  '  C  '  matrix 
do  120  i=l,B 

if  (NF(i)  .LT.  N)  then 

C(NF(i),i)=l 
end  if 

if(  NT(i)  .LT.  N)  then 
C(NT(i),i)=-l 
endif 
c     convert  diameter  in  inches  to  feet 
d(i)=d(i)/12.0 
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area(i)=pi*(d(i)**2)/4 

length=ell(i) 
c    60  diameters  accounts  for  in/out  head  loss 

addl=RLCF*d(i) 

ell(i)=length+  addl 
C   CALCULATE  THE  LAMINAR  RESISTANCES  FOR  EACH  BRANCH 

RLAM(i)=  32.0*  mu*  ell(i)/(AREA(I)  *  (d(i))**2) 

YLAM(I)=1/RLAM(I) 
120  CONTINUE 
C  WRITE  "C"  MATRIX  TO  SCREEN  AND  RECORD  FOR  CHECK 

WRITE(7,123) 
C      WRITE(*,123) 

123   FORMATC  '  THE  C  MATRIX  BELOW,  TRANSPOSED   '  ) 
DO  126  1=1, B 
C      WRITE(*,125)  (C(J,I),J=1,N1) 
WRITEC7.125)  (C(J,I),J=1,N1) 

125  F0RMAT(8(1X,F8.4)) 

126  CONTINUE 

C  WRITE  THE  LAMINAR  RESISTANCE  AND  ADMITTANCE  VECTORS 

WRITE (7,*)  'THE  LAMINAR  RESISTANCE  AND  ADMITTANCE  VECTORS' 
C      WRITE (*,*)  'THE  LAMINAR  RESISTANCE  AND  ADMITTANCE  VECTORS' 
DO  130  1=1, B 

WRITEC7,*)  RLAM(I),YLAM(I) 
C      WRITEC*,*)  RLAM(I),YLAM(I) 
TURB(I)=1.0 
130   CONTINUE 
C  MATRIX  MULTIPLICATION  OF  MATRICES 

C 
C 

DO  140  J=1,B 

DO  140  K=1,B 

Y(J,K)=0 

IF  (J  .EQ.  K)  THEN 

Y(J,J)  =  YLAM(J) 
ENDIF 
140   CONTINUE 

C  MATRIX  MANIPULATION  SECTION 
c  Nl  rows  of  C 
c  b  columns  of  C 
c  CT  transpose  of  C 

CALL  TRANSP(C,CT,N1,B) 
C    write(*,*)  '  write  transpose  of  matrix  C  ' 
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writeC7,*)  '  write  transpose  of  matrix  C 
call  Warray(CT,B,Nl) 

writeC*,*)  '  write  matrix  Y' 
write (7,*)  '  write  matrix  Y' 
call  Warray(Y.B.B) 


C  PROPOSED  ENTRY  POINT  FOR  ITERATIVE  SCHEME 

250    call  matmul(Cy,C,Y,Nl,B,B) 

C     writeC*,*)  '  write  PRODUCT  C*Y   N1XB  ' 

write (7,*)  '  write  PRODUCT  C*Y   N1XB  ' 

call  WarrayCCy.Nl.B) 

CALL  MATMULCYN,Cy,Ct,Nl,B,Nl) 
c  node  flow  source  vector  CQs 

call  matmuK  Is,C,QS,Nl,B,l) 
c   commented  out  artifact  of  old  scheme 
DO  200  1=1, Nl 

ISCI.D— ISCI.D 

ISI(I,1)=IS(I,1) 
200   CONTINUE 

DO  295  1=1, Nl 
DO  295  J=1,N1 

YNI(I,J)=YN(I,J) 

295  CONTINUE 

C  SECTION  ALSO  WRITES  THE  ARRAY  PRIOR  TO  BEING  cHOLESKY  PROCESSED 

WRITEC7,*)  'YNI  MATRIX  PRIOR  TO  CHOLESKY   ' 
C      WRITEC*,*)  *YNI  MATRIX  PRIOR  TO  CHOLESKY   ' 

call  Warray(YNI,Nl,Nl) 

WRITEC7,*)  'ISI  MATRIX  PRIOR  TO  CHOLESKY   ' 
C      WRITEC*,*)  'ISI  MATRIX  PRIOR  TO  CHOLESKY   ' 

call  Warray(ISI,Nl,l) 

call  CH0LESKY(YNI,ISI,N1) 

WRITE (7,*)  'YNI  MATRIX  FOLLOWING  CHOLESKY   ' 
C      WRITEC*,*)  'YNI  MATRIX  FOLLOWING  CHOLESKY   ' 

call  WarrayCYNI.Nl.Nl) 
C  P  NODE  MATRIX  P=YN"-1*IS 
C  H  BRANCH  PRESSURE  CT*P 
C  Q  BRANCH  FLOW  RATE 
C  YH  MATRIX  PRODUCT  Y*H 

DO  296  1=1, Nl 

P(I,1)«XSI(I,1) 

296  CONTINUE 
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CALL  MATMUL(H,CT,P,B,N1,1) 
CALL  MATMUL(YH,Y,H,B,B,D 
DO  300  1=1, B 
Q(I,1)=YH(I,1)+QS(I,1) 
C   DIVIDE  PRESSURES  IN  LB/FT2  BY  144  TO  CONVERT  TO  PSI 
300   CONTINUE 

c   REYNOLDS  NUMBER  REGIME  DETERMINATION 
DO  5000  1=1, B 

V(I,1)=ABS(  Q(I,1)/AREA(I)) 
REY(I,l)=rho*V(I,l)*d(I)/mu 
REN=REY(I,1) 
C    VRITE(*,*)  '  REYNOLDS  #',  REY(I.l) 
WRITE(7,*)  '  REYNOLDS  •',  REY(I.l) 
if  (REY(I.l)  .LE.  2100.0)  THEN 
Y(I,I)=  YLAM(I) 
TURB(I)=1.0 
TURBAV=1.0 
ERR(I)=0 
F=64. 0/2100.0 
else 

if  (REY(I.l)  .IE.  REDGE  )  THEN 
X1=REDGE-2100.0 
Yl=Fcros(i)  -  64.0/2100.0 

F=  64.0/2100.0  +  Y1*SIN(PI*  (REN-2100)/(2*X1)) 
else 
c        colebrook-white  formula 

sqf= (64. 0/2100.0)**. 5 
1310         save=sqf 

sqf=1.0/(-2.0*loglO(2.51/(Ren*save)  +  rat(i)/3.7)) 
C  IF  (SAVE  .EQ.  SQF)  GOTO  1400 

error=abs(save-sqf )/sqf 
if  (error  .gt.  .001)  goto  1310 
1400         f=sqf**2 
end  if 

TURBF=F*Ren/64.0 
TURBAV= (TURB (I ) +TURBF) /2 . 0 
Y(I,I)=1/(RLAM(I)*TURBAV) 
ERR ( I ) = ABS (TURB AV-TURB  ( I ) ) /turb ( i ) 
TURB(I)=TURBAV 
ENDIF 
C     WRITE(*,*)  I.'FRICTION  FACTOR',  F.'THE  TURB  FACTOR' .TURBAV 
WRITE(7,*)  I.'FRICTION  FACTOR',  F.'THE  TURB  FACTOR' .TURBAV 
5000  CONTINUE 
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EPS=0 

DO  5100  J=1,B 
C     WRITE(*,*)  '  THE  ERROR  FOR  BRANCH', J,'  IS',ERR(J) 

WRITE(7,*)  »  THE  ERROR  FOR  BRANCH', J,'  IS'.ERR(J) 

IF  (ERR(J)  .GT.  EPS  )  EPS=ERR(J) 
5100  CONTINUE 

IF  (EPS  .GT.  MAXERR)  GOTO  250 
C  C  CCCCCCC  NOTE  ERR  DEFINED  IN  TERMS  OF  TURBF  IS  ALREADY  NORMALIZED 
WRITEO.5900)  BR.NAME(l)  ,NAME(2)  ,NAME(3)  ,NAME(4)  , 

♦  NAME(5),NAME(6),NAME(7),NAME(8) 
5900  FORMAT (A3, 8 ( IX, A7)) 

C 

C    READY  TO  BEGIN  OUTPUT  OF  SOLUTION 

C      WRITE(*,*)  'BRANCH, FLOW(Q) .PF(MPSI) ,PT(MPSI) ,  head,  ERR, 

C    4  R£Y#,  TURB  FACTOR,   FLOW(Q-QS),  ' 

WRITE(7,*)  ' BRANCH, FLOW (Q),   presURE  F,   pres  T,  head,  ERR, 
+  REYi,  TURB  FACTOR,  FLOW(Q-Qs)' 

DO  7000  J=1,B 
c     Sources  calculated  in  terms  of  effective  pressure  diff  (  thevenin) 
C  NOTE  THE  FLOW  VALUES  ARE  IN  GAL/MIN 
C    CONVERT  FLOW  VALUES  TO  GAL/MIN 

C   DIVIDE  PRESSURES  IN  LB/FT~2  BY  144  TO  CONVERT  TO  PSI 
C   ADMITTANCE  Y  CONVERTED  TO  (GAL/MIN) /mPS I 
c  *  1000  to  convert  to  mpsi 

PF=P(NF(J),1)/144. 0*10**3 

PT=P(NT(J),1)/144. 0*10**3 

H(j,l)=H(j,l)/144. 0*10**3 

Y( J, J)=Y( J, J)*450. 0*144. 0/10**3 

QS(J,l)=450.0*qS(J,l) 

Q(J,1)=Q(J,1)*450 
C      QG(J)=Q(J,l)-Qs(J,l) 

C      WRITE(*,6000)J,Q(J,l),PF,PT,H(j,l),ERR(j), 
C    ♦  REY(J,1),TURB(J),Y(J,J) 

WRITE(7,6000)J,Q(J,l),PF,PT,H(j,l),ERR(j), 
+  REY(J,1),TURB(J),Y(J,J) 

WRITE(9,6000)J,Q(J,l),PF,PT,H(j,l),ERR(j), 

♦  REY(J,1),TURB(J),Y(J,J) 

6000  F0RMAT(I3,1X,F7.2,3(1X,F7.1),1X,F7.4,2(1X,F7.1),1X,F7.4) 
7000  CONTINUE 


WRITE(*,*) 
WRITE(*,*) 
WRITE(*,*) 
WRITE(*,*) 
WRITE(*,*) 


SEE  FILE  flowOUT.DAT  FOR  OUTPUT        ' 
SEE  FILE  flowDIA.DAT  FOR  CHECKING  INPUTS 
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WRITE(*,*)  •  ****************************************> 

STOP 

END 

subroutine  cholesky(a.b.n) 
C  SOLVES  THE  PROBLEM  AX=B 
C    X  UNKNOWN  N  TUPLE 
C    A  KNOWN  RANK  2  MATRIX  OF  ORDER  N 
C    B  KNOWN  N  TUPLE  (  UPON  ENTRY  ) 
C    B=X  (  UPON  EXIT) 

real  b(40,l) 

real  a(40,40) 

write (7,*)  '  in  cholesky  routine,  matrix  order  n  follows 

call  Warray(a,n,n) 

call  decomp(a.n) 

b(l,l)=b(l,l)/a(l,l) 

do  10  i=2,n 

d=b(i,l) 

il=i-l 

do  5  1=1, il 
5     d=d-a(l,i)*b(l,l) 
10    b(i,l)=d/a(i,i) 

b(n,l)=b(n,l)/a(n,n) 

nl=n-l 

do  30  1=1, nl 

k=n-l 

kl=k+l 

do  20  j=kl,n 
20    b(k,l)=b(k,l)-a(k,j)*b(j,l) 
30    b(k,l)=b(k,l)/a(k,k) 

return 

end 

subroutine  matmuKc  ,a,b,n,m,l) 
C     SOLVES  THE  PROBLEM  C=A*B 
c     C  UNKNOWN  MATRIX  lxn 
c     A  KNOWN  MATRIX  MzN 
c     b  KNOWN  MATRIX  lXm 

real  a(40,40),  b(40,40),  c(40,40) 
C   n     #  columns  of  a  and  c 
cm     #  rows  of  a  and  columns  of  b 
c   1     i  rows  of  b  and  c 

do  25  i=l,n 

do  25  j=l,l 
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c(i,j)=0 

do  25  k=l,m 
25  c(i,j)=  c(i,j)+  a(i,k)*b(k,j) 

return 
end 

subroutine  decomp(a.n) 
real  a(40,40) 
i=l 

write (7,*)  '  now  in  decomp  routine   ' 
write(7,*)  'the  order  of  matrix  a  is',  n 
do  410  ii=l,n 

write(7,35)  (a(ii, jj) ,jj=l,n) 
410   continue 

if  (  a(l,l)  )  1,  1,  3 

1  write(7,2) 
write(7,*)  i,a(i,i) 

2  format (  '  zero  or  negative  radicand  ') 
goto  200 

3  a(l,l)=sqrt(a(l,l)) 
do  10  j=2,n 

10   a(l,j)=a(l,j)/a(l,D 
do  40  i=2,n 
il=i-l 
d=a(i,i) 
do  20  1=1, il 
dhold=d 

20  d=d-a(l,i)*a(l,i) 

if  (  d  .eq.  0)  then 

WRITEC7,*)  '  d  is  zero   ' 
write(7,*)  i,a(l,i),  dhold 

endif 

write(7,*)  '  i,  matrix(i.i)  follow  ' 

write(7,*)  i,  a(i,i) 

if  (a(i,i))  1,1,21 

21  a(i,i)=sqrt(d) 

if  (  i  .eq.  n)  goto  45 
i2=i+l 
do  40  j=i2,n 
d=a(i,j) 
do  30  1=1, il 
30    d«d-a(l,i)*a(l,j) 
40    a(i,j)=d/a(i,i) 
45    do  50  i=2,n 
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c     zero  out  lower  part  of  matrix 

il=i-l 

do  50  j=l,il 
50    a(i,j)=0 

write (7,*)  'made  it  through  decomp  ,  matrix  follows 

do  400  ii=l,n 

write(7,35)  (a(ii, j j) ,jj=l,n) 
35    format(8(2x,f8.4)) 
400   continue 
200   return 

end 

SUBROUTINE  TRANSP(A,AT,NR,NC) 

REAL  A(40,40),  AT(40,40) 

DO  20  J=1,NC 

DO  20  K=1,NR 
20    AT(J,K)=A(K,J) 

RETURN 

END 

subroutine  warray(a,n,m) 
c  writes  array  to  screen  and  to  diagnostic  file 

real  a(40,40) 
c    n  rows  of  a 
c    m  columns  of  a 

do  20  j=l,n 

write(7,30)  (a(j ,i) ,i=l,m) 
C     write(*,30)  (a(j ,i) ,i=l,m) 
20    continue 
30    format  (8(lx,f8.4)) 

return 

end 
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